Libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))
Parameters and
risk
censoredProb <- 0.0002
timeSpan <- 20
timeInterval = 0.01
InitialPopulatoin <- 1000
ContBetaRate_1 <- 0.0002
ContBetaRate_2 <- 0.000001
BinBetaRate_1 <- 0.001
BinBetaRate_2 <- 0.003
BaselineHazard <- ContBetaRate_1
betaRates <- c(BaselineHazard,ContBetaRate_1,ContBetaRate_2,BinBetaRate_1,BinBetaRate_2)
BaseVar <- rep(1,InitialPopulatoin)
ContVar_1 <- runif(InitialPopulatoin)
summary(ContVar_1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000764 0.2350957 0.4777109 0.4771627 0.7061444 0.9996590
ContVar_2 <- rnorm(InitialPopulatoin,50,10)
ContVar_2[ContVar_2 < 1] <- 1
summary(ContVar_2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 20.34 44.05 50.62 50.50 57.16 82.34
BinVar_1 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_1)
## BinVar_1
## 0 1
## 628 372
BinVar_2 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_2)
## BinVar_2
## 0 1
## 661 339
dataFeatures <- as.matrix(cbind(BaseVar,ContVar_1,ContVar_2,BinVar_1,BinVar_2))
hazardRate <- as.numeric(dataFeatures %*% betaRates)
summary(hazardRate)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002367 0.0003613 0.0013161 0.0017349 0.0033288 0.0044498
Getting the events
and time to event
aliveSet <- c(1:InitialPopulatoin)
eventSet <- numeric(InitialPopulatoin)
timetoEvent <- numeric(InitialPopulatoin)
for (time in c(1:(timeSpan/timeInterval)))
{
randProb <- runif(length(aliveSet))
Iscensored <- randProb <= censoredProb
Isevent <- randProb <= (1.0-exp(-hazardRate[aliveSet]))
Iscensored <- Iscensored & !Isevent
eventSet[aliveSet] <- Isevent
timetoEvent[aliveSet] <- time*timeInterval-timeInterval/2
isCensoredOrEvent <- Iscensored | Isevent
aliveSet <- aliveSet[!isCensoredOrEvent]
# cat(length(aliveSet),"(",sum(isCensoredOrEvent),",",sum(Isevent),",",sum(Iscensored),")\n")
}
timetoEvent[aliveSet] <- time*timeInterval + timeInterval/2
summary(timetoEvent)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.005 1.915 5.685 8.504 16.413 20.005
table(eventSet)
## eventSet
## 0 1
## 221 779
pevent <- (1.0-exp(-hazardRate))
summary(pevent)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002366 0.0003612 0.0013152 0.0017323 0.0033233 0.0044399
simulatedDataFrame <- as.data.frame(cbind(status=eventSet,time=timetoEvent,pevent=pevent,dataFeatures))
RRplots()
plotTimeInterval <- 2.0
hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)
rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04623 0.06971 0.23143 0.26261
0.48612 0.58933
table(simulatedDataFrame$status)
0 1 221 779
RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90,0.80),
timetoEvent=simulatedDataFrame$time,
title="Simulation",
ysurvlim=c(0.00,1.0),
riskTimeInterval=plotTimeInterval)






par(op)
By Risk
Categories
obsexp <- RRAnalysisCI$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
779 |
725.2 |
835.7 |
770.5 |
1.011 |
0.941 |
1.08 |
0.759 |
| low |
178 |
152.8 |
206.2 |
167.8 |
1.061 |
0.910 |
1.23 |
0.418 |
| 90% |
29 |
19.4 |
41.6 |
28.1 |
1.032 |
0.691 |
1.48 |
0.850 |
| 80% |
572 |
526.1 |
620.9 |
575.6 |
0.994 |
0.914 |
1.08 |
0.900 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event
Analysis
isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Observed Time",
at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
| 1Q |
3.57 |
3.25 |
| Median |
9.38 |
5.44 |
| 3Q |
13.34 |
11.45 |
Table continues below
| 1Q |
1.13 |
13.4 |
| Median |
2.56 |
15.1 |
| 3Q |
5.71 |
16.6 |
| 1Q |
11.4 |
1.46 |
| Median |
11.5 |
1.52 |
| 3Q |
11.7 |
3.66 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

Risk
Calibration
op <- par(no.readonly = TRUE)
crdata <- cbind(simulatedDataFrame$status,pvalue,simulatedDataFrame$time)
#calprob <- CalibrationProbPoissonRisk(crdata,timeInterval=plotTimeInterval)
calprob <- CalibrationProbPoissonRisk(crdata)
( 20.28598 , 11004.11 , 1.935431 , 779 , 1494.978 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
By Risk
Categories
obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
779 |
725.2 |
835.7 |
772.4 |
1.009 |
0.939 |
1.08 |
0.801 |
| low |
178 |
152.8 |
206.2 |
168.2 |
1.058 |
0.908 |
1.23 |
0.441 |
| 90% |
29 |
19.4 |
41.6 |
28.2 |
1.030 |
0.690 |
1.48 |
0.850 |
| 80% |
572 |
526.1 |
620.9 |
577.0 |
0.991 |
0.912 |
1.08 |
0.851 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event
Analysis
isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Observed Time",
at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
| 1Q |
3.57 |
3.25 |
| Median |
9.38 |
5.44 |
| 3Q |
13.34 |
11.45 |
Table continues below
| 1Q |
1.13 |
13.4 |
| Median |
2.56 |
15.1 |
| 3Q |
5.71 |
16.6 |
| 1Q |
11.4 |
1.46 |
| Median |
11.5 |
1.52 |
| 3Q |
11.7 |
3.66 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)
